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Abstract 

We study the approach to equihbrium of a Bose gas to a superfluid state. We 
point out that dynamic scahng, characteristic of far from equihbrium phase- 
ordering systems, should hold. We stress the importance of a non-dissipative 
Josephson precession term in driving the system to a new universality class. 
A model of coarsening in dimension d = 2, involving a quench between two 
temperatures below the equilibrium superfluid transition temperature (Tc), is 
exactly solved and demonstrates the relevance of the Josephson term. Nu- 
merical results on quenches from above Tc in d = 2,3 provide evidence for the 
scaling picture postulated. 
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The recent observation of Bose condensation in neutral, trapped atomic gases [H and 
excitons in CU2O opens up exciting possibilities on experimental studies of time-dependent 
non-equilibrium phenomena in a heretofore inaccessible regime. In particular, an issue which 
could be experimentally investigated, and which we shall address theoretically in this paper 
is the following — upon quenching a Bose gas to a final temperature (T) below T^, how does 
the condensate density grow with time before attaining its final equilibrium value? A few 
recent papers have addressed just this question, but they have focussed on the early time 
(on the order of a few collision times), non- universal dynamics. However, as has also been 
noted recently in Ref , the interesting experimental questions are instead associated with 
the long-time dynamics involving the coarsening of the Bose condensate order parameter. 
This dynamics is "universal" in a sense that will be clarified below. 

A natural and precise language for describing the evolution of the condensate is offered 
by recent developments in the theory of phase-ordering dynamics in dissipative classical spin 
systems, as reviewed in the article by Bray . In this theory, one considers the evolution of a 
classical spin system after a rapid quench from some high T to a low T in the ordered phase. 
The dynamics is assumed to be purely relaxational, and each spin simply moves along the 
steepest downhill direction in its instantaneous energy landscape. Locally ordered regions 
will appear immediately after the quench, but the orientation of the spins in each region will 
be different. The coarsening process is then one of alignment of neighboring regions, usually 
controlled by the motion and annihilation of defects (domain walls for Ising spins, vortices 
for XY spins etc.). A key step in the theory is the introduction of a single length scale, 
/(t), a monotonically increasing function of the time t, which is about the size of a typical 
ordered domain at time t. Provided l{t) is greater than microscopic length scales, like the 
range of interactions or the lattice spacing, it is believed that the late stage morphology of 
the system is completely characterized by l{t), and is independent of microscopic details, 
i. e. it is universal. This morphology is characterized by various time dependent correlation 
functions which exhibit universal scaling behavior. 

We turn then to the Bose gas. The order parameter in this case is the boson annihilation 
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field (where r is a spatial co-ordinate); the phase of the expectation value of is 

aligned across the system in the equilibrium Bose-condensed state. A key point is that after 
relatively few atomic collisions, once the domain size l{t) is large enough {e.g. larger than 
the de Broglie wavelength), it is permissible 0] to treat ip{r,t) as a classical field which 
obeys Hamilton- Jacobi equations of motion (for a related discussion on the emergence of 
classical dynamics in the equilibrium properties of an antiferromagnet, see Ref 0). It must 
be kept in mind that it is only the equations of motion for the collective order parameter 
which are classical — the very existence of the complex order parameter is due entirely to 
quantum mechanics, and the fact that there is a wavefunction for the condensate. 

An important property of the equations of motion for ip, discussed below, is that they 
are not simply relaxational. Instead, they contain non-dissipative, kinematical "streaming" 
or "Poisson bracket" terms . One such term is responsible for the Josephson precession of 
the phase of ^ at a rate determined by the local chemical potential. A central objective of 
this paper is to understand the consequences of such terms on the phase-ordering theories 
of Ref 1^. We will argue that the Josephson term constitutes a relevant perturbation on 
the dynamics and that the resulting coarsening process belongs to a new universality class. 
Specifically, in the remainder of the paper we will (i) exactly solve a model of a = 2 Bose 
gas always in contact with a reservoir, where the temperature of the reservoir is suddenly 
switched between two temperature below |]10|; {H) present numerical results on the time 
evolution of an isolated Bose gas in d = 2, 3, where the initial state has no superfluid fraction, 
while the final state is superfluid. 

We will begin by considering a solvable coarsening problem which illustrates the possible 
consequences of the Josephson term in a simple setting. A d = 2 Bose gas is superfluid for 
T < TxT, the well known Kosterlitz-Thouless phase transition temperature; consider the 
phase ordering process in which the Bose gas is moved at time t = from contact with a 
reservoir at an initial T = Tj, to a reservoir with a final T = Tf, such that Tf < Ti < Tkt] 



a similar quench was considered in Ref |T^ for the purely dissipative XY model. In the 



long-time limit, all vortices and fiuctuations in the amplitude of ip can be neglected, and we 



may parametrize if) = e**^. The free energy density in the purely dissipative XY model is 
now determined simply by the gradients of the phase ~ (V0)^. In the case of the Bose gas, it 
is also necessary to take the conserved number density into account. Let m be proportional 
to the deviation of the particle density from its mean value; then the free energy density we 
shall work with is 

^=l|rfV[(V0)2 + m2]. (1) 

We have rescaled spatial co-ordinates and m to obtain convenient coefficients in JF. The 
Josephson precession term, whose effects we wish to study, is contained in the Poisson 
bracket 

{'m{r),(f){r')} = go5{r - r'), (2) 

where go is a constant. The method reviewed in Ref now leads to the /mear equations of 
motion HIT 



^ = roV'0 + gom + e, 

ot 

din 

— = Ao V'm + goV^<l) + ( (3) 

where the coefficients Fq, Aq > represent the dissipation arising from the coupling of the 
system to the reservoir. The effects of the reservoir are also contained in the Gaussian 
noise sources 6 and ( with zero mean and (for t > 0) correlations appropriate to T = Ty: 
{e{r, t)e{r', t')) = 2roTf6{r - r')5{t - t'), {({r, t)C{r', t')) = -2XoTfV^6{r - r')6{t - t') and 
{C,{r,t)0{r' ,t')) = {kB = 1). The equations (0) are linear, can be easily integrated, and all 
correlations can be computed exactly. 

Let us first recall the structure of the solutions expected from naive scaling for > 2. 
For our models we expect the domain size l{t) ~ t^/^ where z is a non-equilibrium dynamic 
exponent. We consider the behavior of two correlation functions: {i) The equal-time cor- 
relator G{r,t) = {ilj*{r,t)ilj{0,t)) (the k = component of the spatial Fourier transform of 
G is proportional to the condensate fraction). Under scaling we expect for large r and t, 



G(r,t) ~ r~^f gij- /t^l^) where (7 is a universal scaling function, rjf = for d > 2, while in 
d = 2, r]f (= Tf/2'K for JF^) is the equilibrium exponent of the quasi-long-range order at 
T = Tf. {ii) The unequal time correlation function C{r,t) = {ilj*{r,t)tp {0,0)) for which we 
expect for large r and t, C{r,t) ~ t"^^^ f{r /t^/^) where / is a universal scaling function, and 
A is a new dynamic exponent. 

It turns out that our model JF does not completely obey the simple scaling hypotheses 
as stated above. This becomes clear upon considering the two-time correlation, C, which 
turns out to depend upon two large length scales, li{t) ~ (at)^/^ and l2{t) ~ g^t (with 
a = (Fo + A)/2): it obeys the scaling form C(r,t) ~ t-(3'?»+'?/)/4/(r/(at)i/2, r/(^oi)) (where 
Tji = Ti/2'n). The dependence of these scales on qq suggests that is a relevant perturbation 
with renormalization group eigenvalue 1, in the language of Ref The scaling function / 
was determined to be 

/(xi, X2) = eM-V^ r^{l- My)} cos(y/x2)e"^'/^?]. (4) 
Jo y 

For r ~ hit), we use f{xi,X2 ^ 0) = 1: this shows that the autocorrelation C{0,t) ~ 
t'C^Vr+Vf)/^ in contrast to the resuh in the model of Ref |g C{0,t) ~ On the 

contrary, one could insist on a scaling picture using only the single larger length scale r ~ 
hit), and would then need /(xi — > 00, X2) which equals [l + y^l — x^]^'^^ for X2 < 1 and equals 
X2^" for X2 > 1. It can also be checked that one recovers the initial equal time equilibrium 
result for C(r, t) when r ^ 00 with t large but fixed. We also note that the relevance of go 
is evident in the autocorrelations of m. We find (m(0, t)m(0, 0)) ~ il/t)fiigQ^t/a) where 

/i(r) = A7f\[l - r smye-y'/'^'dy]; (5) 
Jo 

clearly, for g^ = 0, this autocorrelator decays as 1/t for large t, while for nonzero go it 
decays faster as Finally, results on the equal-time ip correlator, G. It has a crossover 
at a time ti ~ ci/gQ with d = |Fo — Ao|; this time is similar to the crossover time in 
(m(0, )f:)m(0, 0)), except that d has replaced a. Both for t <^ ti and for t ^ ti, G obeys a 
scaling form similar to that obtained in the relaxational model of Ref [0 (which has go = 0): 
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G(r,t) ~ r~'^f g{r /(■yty/'^) where g is a. scaling function described in Ref [0]; however, the 
rate 7 = Fq for t ^ ti and 7 = a for t ^ ti. 

While this phase only model- JF is not relevant for studying quenches from above the 
transition temperature (since it neglects the non-linear terms), the exact solution of this 
linear model is quite instructive. It clearly emphasizes the importance of the nondissipative 
Josephson coupling term. In fact as seen above, the presence of this term [qq ^ 0) changes 
the universality class of the system. Thus it is reasonable to expect that even for quenches 
from above the transition temperature, this term would play an important role. In this 
case, after the quench the system has defects (e.g., vortices in 2-d). As time progresses, 
these defects move around and annihilate each other and the system becomes more and 
more ordered. To study this coarsening process that proceeds via the annealing of defects it 
is necessary to study the evolution of both the phase and amplitude of ip. This growth of long 
range order in the system can be studied in two different ways. In one case one considers 
the time evolution of an isolated Bose gas, not in contact with a reservoir. Though the 
dynamics in this case is nondissipative, the system still exhibits irreversible approach to 
equilibrium. In the other case, the Bose gas is in contact with a heat bath. These are the 
analogues of micro-canonical and canonical ensembles in equilibrium statistical mechanics. 
It is reasonable to expect that both descriptions would lead to the same results for universal 
scaling properties. Most previous studies on coarsening have been done in the "canonical" 
ensemble. However in this paper, we use the "microcanonical" approach. To the best of our 
knowledge, this approach has never been used before to study coarsening in any system. As 
we will see below, the dynamics in the "microcanonical" ensemble is completely specified 
by the Hamiltonian of the system with no additional phenomenological parameters. The 
"canonical" dynamics, on the other hand, needs several phenomenological constants as input 
parameters. 

For the isolated Bose gas ("microcanonical" ensemble), an excellent approximation for 
the total energy of an order parameter configuration ip{r^t) is Ti. = J d'^r jVV'P + f IV'I^ 
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where we have rescaled lengths to make the coefficient of the gradient term unity, and m > 
is the two-particle T-matrix at low momentum. The Hamilton- Jacobi equation of motion 
for ^|J follows from the Poisson bracket {tp, tp*} = i 

z^ = [-V' + um^, (6) 

and is the well-known ||12| Gross-Pitaevski (GP) or non-linear Schroedinger equation. We 
can also add a quadratic {ipl"^ term to Ti, and it leads to term linear in ip in the GP equation; 
however this linear term can be eliminated by an innocuous global phase change in ip. The 
GP equation conserves the total number of particles / d'^r\ilj\'^ , the total momentum, and 
7i, and hence there is no global dissipation of energy. Nevertheless, in the thermodynamic 
limit, the GP equation does display irreversible coarsening, as will be abundantly clear from 
our numerical results to be described later: a random initial state with a negligible number 
of particles in the zero momentum (k) state (i.e., short range initial correlations), evolves 
eventually to a state with a condensate fraction equal to that expected at equilibrium in the 
microcanonical ensemble at the total energy of the initial state. 

In the "canonical" approach on the other hand, it is permissible to add dissipative terms 
to the equation of motion ofip. A simple additional damping term to the GP equation leads 
to a model expected to be in the same universality class of the so-called Model A P,^; this 
model is however not acceptable: it violates local conservation of the particle density, and, 
as discussed near (|]), it is necessary ||8|JT3[| to introduce the density fluctuation field, m(r, t); 
the value of \ip{r,t)\'^ is then the contribution to the particle density from low momentum 
states, while m(r, t) represents the density fluctuation from all states; the Poisson bracket in 
this case is {m{r),ilj{r')} = igQilj{r)6{r — r'). This is model F in the language of Ref P| (see 
also 1^). Note that the strength of the crucial precession term in the dynamics is controlled 
by which is an adjustable phenomenological parameter (however, in the Hamiltonian 
dynamics considered earlier, there is no such freedom). Numerical study of coarsening using 
model F could thus be complicated by crossover effects associated with the adjustable value 
of fi'o (fl'o = corresponds to the purely dissipative model-A dynamics, which is clearly in a 



different universality class). 

We therefore restrict our numerical study here to the GP model. All of the numerical 
results obtained so far are consistent with the simplest naive scaling hypotheses described 
earlier, and do not require the introduction of two length scales, as was necessary in the 
linear model above (though we have not yet obtained numerical results on unequal time 
correlations, for which the linear model T clearly displayed two length scales). We will 
present results both m. d = 2 and d = ?>. The d = 2 system allowed us to study larger sizes 
with better finite-size scaling properties. 

We discretized @ on a lattice, and integrated in time using a FFT-based algorithm 
which conserved energy and particle number to a high accuracy. We work in units where 
the lattice spacing is unity and choose the scale of the lattice field to make the density one. 
We set u to be approximately 0.25 so that we are considering a dilute gas. We choose an 
ensemble of initial conditions with a narrow distribution of energy, whose width goes to zero 
in the thermodynamic limit. We assign initial values to the Fourier components 0) as 



follows: '?/'(/;;, 0) = Yno(A;) exp[z</)(A;)] where the (/)(fc)'s are independent random variables 
chosen from a uniform distribution with range [0, 27r] and the function n^ilz) is chosen 
to ensure that initial real-space correlations are short-ranged (corresponding to a "high 
temperature" configuration) while still having low enough energy so that the equilibrium 
state corresponding to this energy is superfiuid. Though the ensemble of initial conditions 
is not strictly the Gibbs distribution at any temperature, it is however expected that the 
precise details of the initial conditions do not matter for the late time universal properties 
as long as the initial correlations are short ranged. 

More specifically we chose no(^) = c/ ((e(A;) -|- /ii)(exp[(e(A;) — ii2)lT\ -|- 1)) where e(A;) is 
the fourier representation of the lattice version of the Laplacian and c sets the overall scale 
of UQ^k). Here one chooses the parameters /ii, fi2 and T to achieve the appropriate trade-off 
between energy and correlation length. Note that this careful choice of initial conditions is 
needed as the GP equation does not have any explicit dissipation and the system evolves in 
phase space on a constant energy surface. 
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We used finite-size scaling to model the results in a finite system of linear dimension L: 
it predicts a scaling form G(r, t) = L"''$G[r/L, t/L^] for the equal time correlation function. 
In d = 3 the exponent r] = 0, while in = 2, it is associated with the final equilibrium state 
and varies continuously with temperature. The structure factor, S{k,t), is obtained by a 
spatial Fourier transform of G{r, t), and the number of particles in the k = mode is clearly 
S{0,t); the latter should satisfy S{0,t) ~ L^-'i^t/L'] in d = 2 and S{0,t) ~ L^^t/L'] 
in d = 3. The scaling function $ goes to a constant for t ^ and the system attains 
equilibrium after a time t ^ L^. 

Results for d = 2 are shown in Fig 0. We performed finite size scaling analysis for 
L = 16, 32 and 64 and found reasonable data collapse with t] ^ 0.27 and z ~ 1.1. The 
value of 7] indicates that we are at a non-zero temperature close to T^t] strictly speaking 
we must have t] < 1/4, but the value of t] is relatively T independent near T^t, and the 
discrepancy is within our numerical errors. The value of z is in sharp contrast to the 
z = 2 (with logarithmic corrections) result obtained by various groups p^ , p!5[| for the purely 
dissipative Model A dynamics ^ (obtained from Model F by setting go = and ignoring 
m) of classical XY spins. Although we have determined the value of z for a quench to a 
particular temperature Ty, we expect that z is the same for all < Tj < T^t- Results for 
d = 3 are shown in Fig |^ for linear sizes L = 16, 32. The data collapse is not as good as 
that in d = 2, but again we obtained a z ~ 1.1. 

Thus our numerical results, both in d = 2 and 3, are consistent with a value of ^ = 1, 
which is also the result suggested by the exact calculation in the phase only model. 

Finally, we close with some physical discussion on reasons for the difference between the 



GP model, and quenches in the purely dissipative Model A [p^p!5[ . The dynamics in the 
GP model proceeds via the annihilation of nearby vortex-antivortex pairs (in = 2) as in 
Model A. However there is an important difference between the two in the details of the 
vortex motion. In Model A, oppositely charged vortices attract each other with a force that 
falls off as the inverse of their separation (apart from logarithmic corrections). Since the 
dynamics is overdamped, this implies l{t) ~ t^^^. In the GP model, on the other hand, 
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the situation is more complex. In addition to vortices, the system also has a propagating 
"spin-wave" mode arising from the streaming terms. A pair of oppositely charged vortices, 
apart from attracting each other, also interacts with the spin-wave background. In addition, 
it experiences a Magnus force which causes the pair to move with uniform velocity in a 
direction perpendicular to the line joining them. These qualitative differences in the nature 
of the defect dynamics change the universality class of the coarsening process. 

In summary, we have presented evidence, both analytical and numerical, that the phase- 
ordering dynamics of the Bose gas belongs to a new universality class. A particular conclu- 
sion of this work is that the condensate density of the Bose gas, following a sudden quench 
from the normal to the superfluid phase in dimensions d >2, will grow at late times as i'^/^. 
We have presented evidence, both analytical and numerical, that z = 1. 

We thank D. Kleppner for stimulating our interest in this problem, and M.P.A. Fisher 
and A. Bhattacharya for useful discussions. This research was supported by NSF Grants 
DMR-92-24290 and DMR-91-20525. 
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FIGURES 




FIG. 1. Numerical results from the simulation of the GP equation in d = 2. The number of 
particles in the zero momentum state is S{0, t) and the figure shows its scaling properties as a 
function of system size, L, and time, t. The inset shows the scaling of the equilibrium equal time 
correlation function, G{r, t oo). The best scaling collapse was obtained in both plots for r] ~ 0.27 
and z ~ 1.1. The scale of all axes (except the values of r/L) are arbitrary. 
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FIG. 2. Numerical results for the GP equation in d = 3. The notation is as in Fig. ^ with 
exponent z ~ 1.15. 
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